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Abstract 

A granular gas may be modeled as a set of hard-spheres undergoing inelastic collisions; its microscopic dynamics 
is thus strongly irreversible. As pointed out in several experimental works bearing on turbulent flows or granular 
materials, the power injected in a dissipative system to sustain a steady-state over an asymptotically large time 
window is a central observable. We describe an analytic approach allowing us to determine the full distribution of 
the power injected in a granular gas within a steady-state resulting from subjecting each particle independently 
either to a random force (stochastic thermostat) or to a deterministic force proportional to its velocity (Gaussian 
thermostat). We provide an analysis of our results in the light of the relevance, for other types of systems, of the 
injected power to fluctuation relations. 

Resume 

Puissance injectee dans un gaz granulaire. Un gaz granulaire est assimile a une assemblee de spheres dures 
inelastiques, ce qui confere a sa dynamique microscopique un caractere fortement irreversible. Ainsi que suggere 
par plusieurs travaux experimentaux, realises sur des fluides turbulents ou granulaires, la puissance injectee dans 
ces systemes pour les maintenir dans un etat stationnaire sur une grande fenetre temporelle est une observable de 
choix. Nous presentons une methode analytique permettant d'acceder a la distribution de la puissance injectee dans 
un gaz granulaire au sein d'un etat stationnaire obtenu en soumettant chaque particule soit a une force aleatoire 
(thermostat stochastique) , soit a une force deterministe proportionnelle a sa vitesse (thermostat gaussien). Nous 
analysons nos resultats a la lumiere du role joue, dans d'autres types de systemes, par la puissance injectee dans 
le contexte des relations de fluctuation. 
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1. Motivations 

One of the lessons of thermodynamics is that in order to characterize the macroscopic properties of systems 
in equilibrium only a few numbers of well-chosen variables are necessary. Intensive state functions, like the free 
energy per unit volume, depending on intensive variables only (like temperature and density), can then be used 
to determine e.g. the system's phase diagram. Statistical mechanics is the theory that precisely allows to bridge 
the microscopies to the interesting macroscopic behavior, which it does by bypassing the details of the dynamical 
rules governing the time evolution of the system at hand. The outcome of its machinery lies in the determination 
of the much sought after state functions. 

Unfortunately, for systems that are in steady-states without being in equilibrium, no such theory exists. As 
long as the dynamics breaks the detailed balance condition, steady-state properties will strongly depend on the 
details of the microscopic evolution rules, and any statistical mechanical approach must incorporate those. Yet, 
a few guiding principles inherited from equilibrium thermodynamics can be saved. We shall illustrate this on the 
example of the granular gas, that we model as an assembly of hard-spheres interacting only through inelastic 
collisions. The irreversible dissipation of energy is compensated by an energy injection mechanism -two such will 
be presented in the sequel- that maintains the gas in a steady-state, with time translation invariant properties. 
Our interest goes to the total energy W{t) injected by the heating mechanism over an asymptotically large time 
window [0, t]. This is a space integrated quantity, as it takes into account the work performed on each individual 
particle over [0,t]. That W is a space extensive quantity means that its properties will not be too sensitive to the 
microscopic details of the interactions. The key-role of the dynamics is taken into account in the time extensivity 
of W , but again with the hope that irrelevant details will be smoothened out through time integration. These 
arguments have already been put forward in a series of experimental and theoretical works bearing on turbulent 
flows, convection experiments, or on granular materials themselves. The belief is that global, i.e. space and time 
integrated, observables, more than local response functions, will allow comparing systems as far apart as turbulent 
flows or granular gases whose dynamics nevertheless share strongly dissipative features, and that the theorist craves 
to unite within a common framework. 

In the present work we shall focus on the probability distribution function (pdf) P(W, t) of W, as well as on 
its related generating function P(X,t) = {e~ xw }. More specifically, our efforts will be devoted to determining the 
corresponding large deviation functions defined by 



n(w) = lim -\nP{W = w t,t), n(X) = lim \ In P(X, t) (1) 

t^oo t t — >oo t 



As has been recognized long ago [1], tt(w), or fi(X), which are related by a Legendre transform tt(w) = max.\{fi(X) + 
Xw}, play the role of intensive dynamical free energies, solely depending on intensive variables, like w(t) = ^j^-, 
the fluctuating time averaged injected work. 

It also turns out that in other classes of dynamical systems, the power injected by the heating mechanism can 
be related to the entropy current flowing into the system, whose large deviation function (ldf) has been shown 
to possess an important symmetry property: this is the celebrated fluctuation theorem [2,3]. While a crucial 
hypothesis underlying its demonstration in the aforementioned works is missing in the granular gas, namely a 
weak form of time reversibility, empirical attempts have been made in the past, at an experimental [4] and a 
theoretical [5,6] level, to see whether such a fluctuation relation for W(t) would hold. There is furthermore no 
straightforward connection between W and entropy production, an otherwise ill-defined concept [7,8] in a granular 
gas. These issues will be further discussed in the final section. For now we start by describing the dynamical 
evolution rules of our granular gas, along with two possible energy injection mechanisms. Then we present a 
kinetic-theory based approach to determining the related large deviation functions. The discussion will cast our 
results within the framework of "fluctuation theorems". 
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2. Stationary state of a granular gas 



2.1. The microscopic dynamics of a granular gas is irreversible 

A standard way of modeling a granular gas is to consider a set of N hard-spheres undergoing inelastic collisions 
in which a fraction (1 — a 2 ) of the relative kinetic energy is dissipated away. The restitution coefficient a lies 
between and 1 (elastic collisions). Two incoming particles with velocities vi and V2 acquire velocities v{ and 
V2 after having collided, with 

Vi = vi - 1 + a (V12 ■ cr)(T, V2 = v 2 + 1 Q (vi 2 • er)cr (2) 

The component of the relative velocity V12 — vi — v-z along the unit vector er joining the centers of the particles 
is reduced by a factor a, v$ 2 ■ & = — avi2 • cr. For any given trajectory in phase space, the time-reversed one is 
not a physical trajectory. One can write an energy balance equation for the total kinetic energy E(t) = 5^ §v 2 , 
which varies according to 

AE = E(t) - E(0) = W(t) - D(t) (3) 

where W(t) is the energy injected by an external driving mechanism, while D(t) > is the energy irreversibly 
dissipated through the inelastic collisions. The average variation rate of D can be estimated as the collision rate 
times the energy dissipated through a collision, 

where I is the mean free path. The mean kinetic energy per particle provides a typical energy scale, also termed 
granular temperature, and it is defined as 

T g = (v?)/d = ft 1 (5) 

In the absence of a heating mechanism, dimensional analysis and the assumption of homogeneity, lead to T g 
decaying with time as T g (t) oc t~ 2 (Haff's law). We now introduce two relevant heating mechanisms. 

2.2. Free cooling and deterministic isokinetic thermostat 

Each particle i, in addition to the inelastic collisions, is subjected to an energy-injecting viscous friction force 
+7Vi. However unphysical it may appear at first sight, this thermostat is actually relevant to the study of the 
homogeneous cooling regime [9,10]. During the homogeneous cooling stage, the typical velocity \fTg, and thus 
the collision frequency, decrease as 1/t. Rescaling time with the collision frequency allows to eliminate time 
dependence and leads to an effective (nonequilibrium) steady-state. The latter rescaling is exactly accounted for 
by an effective viscous friction force that, instead of dissipating energy, pumps it into the system. Within the 
newly rescaled dynamics, collisions occur at a constant rate. The energy provided by the thermostat reads 

w ( t ) = 52 fdr 7 v?(T) (6) 
i J o 

Note that W > for all trajectories in phase space. 



2.3. Heated gas and stochastic thermostat 

A more conventional way of achieving a steady-state is to inject energy by means of independent random forces 
acting on each individual particle: 

^ = Fi + collisions, (Ff(t)F?(t')) = 2V8 lJ 8 af) 8{t - t') (7) 
at J 

This heating mechanism is easier to handle mathematically than more realistic thermostats because it leads to 
a uniform system, by contrast to boundary drives. It should be noted that there are some experimental setups 
that do achieve a uniform heating [11]. Besides, when the experiment resorts to vibrating walls, within a large 
subvolume far enough from the boundaries, energy and particles are roughly uniformly distributed [4], which 
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allows for easier comparison, but this cannot be thought of as resulting from an effective bulk heating mechanism. 
In the present case, the work provided by the external heating reads 

W(t) = J2 /Vf^-v^t) (8) 

i J o 

Given that (W)/t = 2dF, the typical energy scale is set to T g = (jt^jt^Y^ ■ 



3. Large deviation function for the injected power 

Let W(t) denote the total power injected into the granular gas by the heating mechanism over the time interval 
[0, t]. We begin with introducing a phase space density p(F, W, t) that counts the number of systems in state Y that 
have accumulated over the time window [0,t] a total work W(t) = W. A generalized Liouville equation can be 
written for p, in which the Liouville operator Cw can be split into a W conserving part C co \\ and one accounting 
for changes in W under the effect of the external injection mechanism Anj: 

dt P = c w (r, w) P = £ inj (r, w) P + c con {r) P (9) 

It is a convenient detour to go first to the Laplace transform of p, p(F,X,t) = J dWe~ p(F, W, t), then rewrite 
the Liouville equation in terms of p, 

dtp = C W (T, \)p = £ inj (r, X)p + £ C oii(r)/5 (10) 

The largest eigenvalue p(X) of Cw (A) thus governs the asymptotic behavior of p, 

p(r,\,t)~C(\)e» Wt p(T,\) (ii) 

where p(T, A) is the (right) eigenvector of Cw(T, A) associated to p(X), and C(A) is the projection of the initial 
state on this eigenvector. We have chosen to normalize p to unity, which, given that it has a definite sign, endows 
it with the meaning of a probability distribution (note however that (10) is a not a Liouville equation since it 
does not conserve probability). We shall provide further insight on this later. In order to determine p and p, we 
project the eigenvalue equation onto the one particle subspace. We shall focus on this isokinetic thermostat for 
which the results presented here are new. We arrive at 

M (A)/ (1) ( Vl , A) = - 7 9 V1 ■ (vi/ (1) (vi,A)) - A 7 ^/ (1) (v 1; A) + coll . (12) 

where in the rhs of (12) the loose notation "coll" is a shorthand for the full collision operator acting on the two 
point function . In order to extract physical information from this equation, we resort to the molecular chaos 
hypothesis, which turns the two-body interaction term into 

coll. = i / dv 2 da(v 12 .a) f^/ (1) (vr,A)/ (1) (vr,A)-/ (1) ( Vl ,A)/ (1) (v 2 ,A)") (13) 

Not surprisingly, one recovers the steady-state velocity pdf equation at A = 0. It is interesting that the process 
encoded in (12) can be read off as the original granular gas dynamics in which additional particles are created 
(or destroyed, according to whether A < or > 0) at a velocity dependent rate At; 2 . The eigenvalue p(X) is then 
interpreted as the population growth rate. This remark was numerically exploited, for somewhat different systems, 
by Giardina, Kurchan and Peliti [12]. The splitting of as a product of independent one particle distributions 
is indeed a molecular chaos hypothesis for the system with the non particle conserving fictitious dynamics. The 
Boltzmann equation toolbox offers many ways to arrive at an expression for /i(A). 

Let us start with A ~ 0, that is in a regime for which W lies in the vicinity of (W). Strictly at A = 0, that is 
when (v, 0) yields the steady-state velocity pdf, Sonine expansions have successfully been used to characterize 
deviations from the Maxwell distribution [13]. These expansions work all the better as space dimension d is high, 
given that the phase space contraction due to the inelastic collisions occurs only along one -among d- space 
direction. The coefficients of the expansion turn out to be functions of the reduced variable (1 — a 2 )d~ 1 , thus 
making explicit that a large d expansion is equivalent to a quasi-elastic limit, hence the success of expanding 
around a Gaussian. What cannot be guessed, however, is that in practice d = 2 or d — 3 already turn out high 
enough dimensions for the Sonine expansions to provide quantitatively reliable results. This provides the necessary 
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— Gaussian Thermostat 

- - Stochastic Thermostat 




Figure 1. Plot of /x(A). 



Figure 2. Plot of -k{w). 



motivation for attempting a Sonine expansion at A ~ 0. However, further analysis, in the same vein as that carried 
out by Van Noije and Ernst [13], shows that in the large velocity limit (and at A > 0), the one point function 
/^(v, A) must be a Gaussian. The simplest approximation is thus to project /^(v, A) onto a Gaussian with a 
A-dependent variance. From (12) one can easily see that 



HW/N = - 7 A<v 2 ) A 



(14) 



where (. . .)\ denotes an average performed with respect to the weight given by f^K This equation relates the 
variance of j' 1 ' to [i in a particularly simple way. Thus, working in terms of rescaled quantities, with the granular 
temperature precisely given by 



and performing the following replacements 



fi:= n/iNi), \~XT g d/2 



one has, for the dimensionless quantities: 



The asymptotic behavior is given by: 



f (A) = 



1 + 2A-V1+4A 
2A 



M (A) A ^~ 1/4 1 - 2Vl +4A + O ( A + 



„(A) A ~°° -1 - 



7T< 



(15) 
(16) 

(17) 
(18) 



'9 2A 



Taking (14) into account one sees that the typical temperature scale as given by Z' 1 ' is given by T(A) ~ T g 

That T(A — > +oo) — -> means that values of W much lower than the average (W) are produced during trajectories 
over which typical velocities are small, which is of course no surprise. By contrast, given that T(A — > —1/4) = 2T g , 
the larger than average W trajectories arise from realizations in which the effective kinetic temperature is at a 
value twice the stationary state temperature T g . 



It is possible to express the large deviation function of w = W/t by a simple Legendre transform, 



n(w) 



1/3 



+ 00 



2/3^ 



ty(w) 



(19) 



but the validity of the w — > + expression may be hindered by subtle effects that we shall discuss in the next 
section. 

We refer the reader to [14,15] for details about the stochastic thermostat for which we obtain, in terms of 
dimensionless variables (for which (w) = 1), 



T 
dFN 



A := AT„ 



(20) 



the following expression 

M(A) = -A+^^A 1! (21) 

where y / T(X) sets the typical velocity scale for trajectories characterized by A. At large values of A, corresponding 
to values of W small with respect to (W), we expect that T(A) <!C T g and indeed ~ j as A — > +00. This can 
be further refined to obtain the behavior of /i(A) at A — > +00, 

/i(A) A ^°° -A 1 / 4 (22) 
Besides, /i(A) possesses a cut in the A plane at A c = — 3/2 8 ^ 3 , such that, as A — > A+, 

A*(A) = 3/2 2/3 - 3 3/2 2 1/6 v/A - X + C(A - A c ) (23) 

The presence of this cut is responsible for the exponential decay of the pdf of W at large values of W. However, 
the non-analytic behavior at A — > +00 leads to a non-analytic behavior as w = W/t — + + , 

tt(to -> + ) ~ -w~ 1/3 , n(w -> 1) ~ -(to - l) 2 /2 , ?r(w -> 00) ~ -w (24) 



4. Dissipated energy 

4.1. Cumulant generating function 

So far our interest has been focused on the distribution of the energy injected into the system to keep it in a 
nonequilibrium steady state. As already mentioned in section 2, the dissipated energy is intimately related to the 
injected power trough Eq (3) , which expresses the energy balance of the system, and which we rewrite here for a 
practical purpose: 

AE(t) = W(t) - D(t) (25) 

Since we are interested in the large time behavior of both W(t) and D(t), it is worthy to note that while W(t) and 
D(t) are of order t, the boundary term AE(t) is of order one. This simple remark would intuitively, and naively, 
lead to the conclusion that at large time W(t) and D(t) are distributed in the same way, and hence share the 
same large deviation function. Before entering the details of the reasons for which W and D do not necessarily 
share the same large deviation function, we will first show that they actually share the same cumulant generating 
function. 

The dissipated energy D(t) increases at each collision by an amount equal to the difference of the energy of the 
colliding pair just before (at t = t~) and immediately after (t = t + ) the collision: 



D{t + ) = D{t-) + ^—^{v l3 -crf (2(>) 



If we now define the joint probability Pd(T, D, t) of having a microscopic configuration V and a given value of D 
at time t, it is possible to write a Liouville equation for po'- 

^Pd(T, D, t) = C D (T, D)p D (T, D, t) = £inj(r)p D + £ diss (T, D)p D (27) 
where d n j describes some energy injection mechanism. As for the collision operator, its explicit expression reads 

Aiiss = o-^ 1 / dcr(vij ■ a) ( -^5(r 43 - a)b** - <5(r 43 + er) j (28) 

where a is the particle diameter and where 

b**g{vi,. . . ,Vi,. . . ,Vj, ...,v N ,D) =5^vi,...,v**,---,v**,...,Viv,D- ^ J_ _ 1^ (29) 

We now translate (27) in terms of the Laplace transform po(T, X,t) = / dDc~ XD pd(F, D, t), which evolves 
according to 

^p D (T,X,t) =£inj(r) / 3 D +£ d iss(r,A)pi3 (30) 
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with 

£ diss (r,A) = a d - 1 ^ f d*(v tJ ■ *)( ±8(r tJ - a^-*^ 11 ^^- 1 ) b** - 5^ + *)) (31) 

i<3 ;J Vir <T>0 \ a J 

and 

b**g (vi, . . . , Vi, . . . , Vj, . . . , vjv) = 3 (vi, . . . , v**, . . . , v", . . . , vjv) (32) 
A straightforward calculation shows that 

r /T^ \\ -\AE r / 1 i r,\ +AAE -AAE ^. ; r \ +AAB 

-Cdiss(r, A) = c £ diss (r,0)e T = e £ d i ss (r)e T (33) 

Assuming, for definiteness, a deterministic energy injection mechanism in which each particle is subjected to a 
force Fi, we must have 

£ inj (r, A) = - A ^2 Fi • vi - 53 9vt ■ ( F < ) (34) 

i i 

Hence, conversely to (33), we may verify that, 

£ inj (r,A) =e+ AAB £ inj (r,0)e- AAB =e+ AAB £ inj (r)e- AAB (35) 

Using that Cw(T, A) = £inj(r, A)+£diss(I\ 0) and that £d(L, A) = £i n j(r, 0)+£diss(r, A) we arrive at the following 
key identity 

Cw(T,X) =e- XAE C D (r,\)e XAE (36) 

That £d(T, A) and £w(T, A) are related through the similarity transformation (36) establishes that these operators 
have exactly the same eigenvalues, and that the corresponding eigenvectors differ only by a factor e AAB . This 
allows us to conclude that W and D have the same cumulant generating function. Thus, if furthermore the 
largest eigenvalue of jCd(X) is related to the large deviation function for the dissipated energy through a Legendre 
transform, the above analysis shows that the large deviation function of the injected and dissipated energy are 
indeed equal. 



4.2. Influence of time boundary contributions 



It has already been stressed [16] that it is possible that IJ,(X) and n(w) may not be related by a Legendre 
transform. As already discussed in the previous subsection, one expects, in writing the energy balance equation 
(25), that, given that AE is not extensive in time, W and D should have the same large deviation functions 
(expressed in terms of w = W/t or 5 = D/t). Since they share the same cumulant generating function /x(A), this 
also means that within a finite interval around their average value, w and d do share the same large deviation 
function. Nevertheless, beyond a finite value of w or 5, their ldf may become distinct; this may occur [16,17,18,19] 
if AE is distributed exponentially or slower than exponentially. The technical reason of this phenomenon is simply 
due to a problem of Laplace transform inversion. In fact, wc assumed that the ldf it is related to the eigenvalue 
\x trough a Legendre transform. This result is obtained carrying out the Laplace transform inversion trough a 
saddle point expansion (in the t — > oo limit). In practice this last step is valid only for values of A for which it 
is possible to define a path (in the A-complex plane) including a straight line parallel to the A-imaginary axis. 
The problem arises hence when there are some cuts which are not included in the expression of /ii(A), but in 
some other subleading term. Unfortunately, as far as we know, there is no general method to know a priori if and 
where the two large deviation functions may differ. The only argument in this direction involves the probability 
distribution of the boundary term: as already mentioned above, if its tails decrease exponentially, or slower, it 
is possible that such problem may occur. This is the scenario in the cases at hand. One can show that AE has 
tails decreasing as exp(-N 1/2 (f3 g AE) 1/2 ) and exp(-A 1/4 (/3 9 AE) 3/4 ) for the Gaussian thermostat and for the 
stochastic thermostat, respectively (an analysis of the typical values of AE is presented in [20]). Note also that 
these tails decrease exponentially in the number of degrees of freedom, which leaves room, in principle, even in 
the thermodynamic limit, for W and D having ldf differing significantly beyond a finite threshold value. 



5. A critical discussion of numerical results 



In this section we would like to present simulations of the granular gas heated with the stochastic thermostat. 
The total work W(t) provided by the random forces over [0, t] is of course, on average, a positive quantity, 
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Figure 3. Plot of tz(w, t) — tt( — w, t) for t — 1, 2 and 3 mean 
free times (mft). The inset show the probability density func- 
tion of w(t) for the same times. 



Figure 4. Plot of the third cumulant of W divided by time. 
The asymptotic regime is reached after fifty mft. 



but there exist phase space trajectories which will yield a negative W (this cannot occur for the deterministic 
thermostat where W{t) > 0). It is therefore tempting to plot the quantity ir(w) — tv(~w) as a function of w. It is 
well-known [2,3] that for a well-defined class of thermostatted systems the power W injected by the thermostat 
verifies the celebrated fluctuation theorem ir(w) — 7r(— w) — (5w, where /3 is an appropriate inverse energy scale. 
However, one of the key hypotheses underlying the latter fluctuation relation is that the dynamics features a weak 
form of time reversibility, for which each trajectory possesses an time reversed partner (however unlikely). In our 
granular gas, dissipative collisions hamper time reversed trajectories to be physically acceptable trajectories at 
all. It should therefore come as no surprise that no specific fluctuation relation emerges in our case. Of course, this 
is fully confirmed by the explicit calculations presented in the previous section, that show no specific symmetry 
property of the injected power ldf. In spite of this, it may be instructive to plot w i— » j In ^ttt^A a t large times (in 
the infinite time limit, this is exactly 7r(«;) — n(—w)), as has been done in a vibrated granular gas experiment [4]. 
This is shown in figure 3, and the result is rather intriguing: one actually observes a straight line with slope (3 g \ 
This deserves to be explained, given that our analytic results establish the absence of such a relation, even for 
values of w not too far from its average. The first remark is that a numerical simulation is always carried out at 
finite times, and thus one is measuring Tv(w,t) — | In P(wt,t) rather than its t — > oo limit 7r(w). Hence the first 
question is: has the simulation reached the infinite time regime? The answer to that question is tricky. Though 
n(w,t) notably deviates from a quadratic form, which would correspond to P(W, t) being Gaussian, this is no 
proof that the asymptotic regime has been reached. For that matter it is instructive to investigate the behavior 
of the third cumulant (W )/t as a function of time. This is presented in figure 4. The conclusion is that the third 
cumulant reaches its asymptotic value, consistent with the analytic expression, at times a few tens as large as 
those for which it is possible to measure negative values of w — W/t, as presented in figure 3. A similar plot of the 
fourth cumulant would signal that the asymptotic regime has not been reached over the chosen time window. What 
one actually observes is simply the leftover of a short time quadratic behavior for n(w). This is also consistent 
with the quadratic approximation for p(A), which, if taken for granted as the whole function, would indeed imply 
tt(w) — tv(—w) — f3 g w (after restoring the appropriate physical scales as in (20)). The lesson to be drawn from 
this discussion is that it seems to be an optimistic endeavor to investigate solely on numerical or experimental 
grounds the validity of the fluctuation relation. Without criteria for the time scale at which the asymptotic regime 
is entered, the analysis is at risk of remaining confined within short time effects [4,21]. 



6. Conclusion 



We have presented an analytic calculation for the large deviation function of an ./V-body dynamical system with 
strong dissipative interactions. By contrast to existing approaches, based on a stochastic modeling, the noise source 
is completely contained within the microscopic formulation, which avoids the arbitrariness inherent to choosing 
a particular type of Markov dynamics. Our calculation required that we extend the standard methods of kinetic 
theory, so as to grasp the infinite hierarchy of correlation functions encoded in a temporal large deviation function. 
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We have shown that theoretical input was required to properly analyze numerical data when large deviations are 
measured. Our study also allows us to infer that the power injected into a granular gas to maintain a steady-state 
has both generic features and ones that are sensitive to the details of the mechanisms at work. At large injected 
powers, the ldf ir(w) behaves linearly with w = W/t, meaning that the power distribution decays exponentially. A 
less robust feature, however, is the w —> behavior, which proves extremely sensitive to the details of the heating 
mechanism. The present work calls for further investigation in other dissipative systems, such as turbulent flows. 
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